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ABSTRACT 

We investigate protostellar collapse of molecular cloud cores by numerical simulations, taking into 
account turbulence and magnetic fields. By using the adaptive mesh refinement technique, the collapse 
is followed over a wide dynamic range from the scale of a turbulent cloud core to that of the first 
core. The cloud core is lumpy in the low density region owing to the turbulence, while it has a 
smooth density distribution in the dense region produced by the collapse. The shape of the dense 
region depends mainly on the mass of the cloud core; a massive cloud core tends to be prolate while 
a less massive cloud core tends to be oblate. In both cases, anisotropy of the dense region increases 
during the isothermal collapse (n < 10^^ cm~^). The minor axis of the dense region is always oriented 
parallel to the local magnetic field. All the models eventually yield spherical first cores (n > 10^^ cm~^) 
supported mainly by the thermal pressure. Most of turbulent cloud cores exhibit protostellar outflows 
around the first cores. These outflows are classified into two types, bipolar and spiral flows, according 
to the morphology of the associated magnetic field. Bipolar flow often appears in the less massive 
cloud core. The rotation axis of the first core is oriented parallel to the local magnetic field for bipolar 
flow, while the orientation of the rotation axis from the global magnetic field depends on the magnetic 
field strength. In spiral flow, the rotation axis is not aligned with the local magnetic field. 

Subject headings: ISM: jets and outflows — ISM: clouds — ISM: magnetic fields — magnetohydrody- 
namics — stars: formation — turbulence 



1. INTRODUCTION 

Magnetic fields and interstellar turbulence are believed 
to play important roles in the gravitational collapse of 
molecular cloud cores. The measured magnetic fields of 
molecular clouds and molecular cloud cores are strong 
and the magn etic energy is approximately equal to the 
kinetic energy (|Crutcher| [r999) . Magnetic fields therefore 
have the potential to control the gravitational collapse 
of cloud cores. Molecular clouds exhibit supersonic line 
widths, which are interpre ted as supersonic turbulence 
(|Zuckerman fc Evanslll973 ). Such supersonic turbulence 
seems to be common in the wid e range from mole cular 
clouds to molecular cloud cores (Langer etJOHMSl). 

One of the important properties of molecular cloud 
cores is their shape, which is likely related to mag- 
netic fields and turbul ence (see, e.g., the review by 
iMcKee fc Ostrikerll2QQ7[ ). Some studies of the shape of 
cloud co res have suggest ed that th ey tend to be pro- 
late (Mv ers et all 1991; R vdenlll996f ). The origin of the 
shape is however still unknown. Protostellar outflow is 
also an important feature related to magnetic fields. Re- 
cent high-resolution observations of submillimeter polar- 
ization succeeded in resolving the magneti c field around 
young stars to ^ 1Q^~^AU (e.g., Henni ng et al. 
iWolf et al.l[200a iVallee et al.l r2003£ iGirart et al.l 



relation between outflow and magnetic field. 

There have been very few theoretical studies of col- 
lapse of magnetized turbulent cloud cores in protostars, 
despite the importance of turbulence and magnetic fields. 
Although self-gravitational turbulent si mulations have 
been performed by many researchers (e .g., iG ammie et al.l 
[20031 : iLi et an[200l lOffner et al.l[2008[ [Bate 20091), most 
investigated large-scale turbulence and focused on cloud 
core formation. A notable exception is the work of 
[Offner et al.[ i 



][^06), 



revealing the structure of the magnetic fields on such 
scales. The outflows often tend to be aligned with the 
m agnetic fleld, whil e some are oriented perpendicular to 
it (|Wolf et al.[l2003[ ). Thus, there is as yet no clear cor- 

Fmatsu@hosei . ac . j pi 

^ Faculty of Humanity and Environment, Hosei University, Fu- 
jimi, Chiyoda-ku, Tokyo 102-8160, Japan 

^ Center for Frontier Science, Chiba University, 1-33, Yayoi- 
cho, Inage-ku, Chiba 263-8522, Japan 



(|2008[ ). who performed high-resolution sim- 
ulations using the adaptive mesh reflnement (AMR) 
technique to study protostellar collapse. However, these 
simulations did not in clude the effects of magnetic flelds. 
[Goodwin et al.[ (|2004f ) also performed high-resolution in- 
vestigations on protostellar collapse of turbulent cores, 
but they also did not take magnetic flelds into account. 

On the other hand, simulations carrie d out to date 
on the formation of the Larson flrst core (|Larson[[T969[ ) 
have taken account of rotation and magnetic flelds but 
not turbulence , both for a ligned rotators ([To misaka[ 
| 2002[: Machida eTal] [20041: [Baneriee fc Pudrit z 2006; 
'Commercon et al. 2010; Tomida et al. 2010) and inclined 
rotators (Mats umo to & To misak a 2004; Machida e t al.[ 
[2006[ : iHennebelle fc Ciardil 120091 ). In these studies, the 
rotation speed and rotation axis are explicitly assumed 
as initial conditions , and the origin of the rotation has not 
been addressed. Burkert fc Bodenh eimerl (|2000|) showed 
that rotation derived from observations could be repro- 
duced by assuming the presence of turbulence with a 
power spectrum of P{k) oc k~'^ with n = 3 — 4. This 
suggests that the rotation of cloud cores originates in 
turbulence. Consequently, taking account of turbulence, 
we can incorporate rotation in the simulations of proto- 
stellar collapse. The present work therefore investigates 
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protostellar collapse of cloud cores, considering both the 
turbulence and magnetic field. The AMR technique is 
adopted in order to resolve the wide dynamic range from 
the cloud core to the first core. 

This paper is organized as follows. In ^ and ^ the 
model and simulation methods are presented. The re- 
sults of the simulations are shown in ^ and they are 
discussed in ^ Finally, some conclusions are given in 

m 

2. MODELS OF CLOUD CORES 

As an initial model of a molecular cloud core, we con- 
sider a turbulent, spherical cloud threaded by a uniform 
magnetic field. The cloud is confined by a uniform am- 
bient gas. 

As a template for a molecular cloud core, we consider 
the d ensity profile of the critical Bonner-Ebert sphere 
(jBonnor 1956; Eb^t 1955). When ^be(0 denotes the 
non-dimensional density profile of the critical Bonnor- 
Ebert sphere (see lChandrasekharlll939f ). the initial den- 
sity distribution is given by 

^(^\ _ / PoQBE{r/a) for r < Rc ... 
^^"^^ ~ \ PoQBEiRc/a) for r > ^ 

and 



where r, G, Cg, and po denote the radius, gravitational 
constant, isothermal sound speed, and initial central den- 
sity, respectively. The gas temperature is assumed to be 
10 K {cs = 0.19 km s~^) . The initial central density is 
set equal to po = 10~^^gcm~^, which corresponds to a 
number density of no = 2.61 x lO^cm"^ for an assumed 
mean molecular weight of 2.3. The non-dimensional pa- 
rameter / denotes the density enhancement factor. The 
critical Bonnor-Ebert sphere is obtained when / = 1. 
An increase in density by a factor / is equivalent to an 
enlargement of the spatial scale by a factor /^/^ for a 
given central density. The radius of the cloud is de- 
fined by Rc = 6.45a = 0.137/^/^ pc, where the numer- 
ical factor 6.45 comes from the non-dimensional radius 
of the critical Bonnor-Ebert sphere. The density con- 
trast of the initial cloud is p{0)/p{Rc) = 14.0. The 
initial freefall timescale at the center of the cloud is 
thus tff = (37r/32Gpo)^/^ = 2.10 x 10^ yr. The mass 
of the cloud core (r < Rc) is M = 2.Slf^'^MQ. The 
spherical cloud described above is located at the center 
of the computational domain of x^y^z G [—2Rc^2Rc] x 
[-2Rc,2Rc] X [-2Rc,2Rc]. 

Turbulence is given as the initial velocity field, and it is 
not driven in course of the simulations, because prestel- 
lar cloud cores are considered here and no driving source 
exists therein. Therefore, free decaying turbulence is con- 
sidered. The initial velocity field is incompressible with 
a power spectrum , of P( k) cx generated according 
to iDubinski et al.l ()1995f ). where k is the wavenumber. 
This power spectrum results in a velocity dispersion of 
(7(A) oc A-^/^, in agreement with the Larson scaling re- 
lations ()Larson| [T981). Note that the scaling relation 
is applied only at the initial condition, and the velocity 
dispersion changes during decay of the turbulence and 
collapse of the cloud cores. The models are constructed 



by changing the mean Mach number of the initial veloc- 
ity field in the range = — 3 with a common template 
of the initial velocity field. 

The initial magnetic field is uniform in the z-direction. 
The field strength is given by Bz = Q^^cr, where 
a denotes the non-dimensional fiux-to-mass ratio, and 
Bcr denotes the critical field strength given by Bcr = 
27rGy ^E (jNakano fc Nakamural 119781 : iTomisaka et"all 
119881 ). The central column density H is calculated by 

E = pdz = 5.38po<^5 where the integral is performed 
along a line passing through the center of the cloud core. 
In this paper, we examine the magnetically supercriti- 
cal core {a < 1). The initial field strength is estimated 
as Bz = 57. 2a/^/^ pG by using the model parameters a 
and /. Note that the model parameter a is the inverse 
of the dimensionless mass-to-flux ratio p {a = 1/p). 

We construct 10 models by changing the three param- 
eters (a, A1, /). Using these parameters, the energies in- 
side the cloud core (r < Rc) are calculated according to 
Appendix[Al The ratio of the thermal energy to the grav- 
itational energy is expressed as £^th/|£^grav| = 0.836/~^. 
When / = 1.86, 3.0, and 6.0, we obtain £^th/|£^grav| = 
0.50, 0.28, and 0.14, respectively. The ratio of the ki- 
netic energy to the gravitational energy is expressed as 
£^kin/|£^grav| = 0.394/~^7W^. For example, the models 
with {MJ) = (1.0,1.68), (3.0,1.68), and (3.0,6.0) ex- 
hibit £;kin/|£^grav| = 0.234, 2.11 and 0.591, respectively. 
Note that the numerical factor 0.394 depends on the seed 
of the random initial velocity field. Finally, the ratio of 
the magnetic energy to the gravitational energy is ex- 
pressed as £^mag/|^grav| = 11.5a^. For example, the 
models with a = 0.1 and 0.25 have £^mag/|^grav| =0.115 
and 0.718, respectively. 

The dynamical evolution of the cloud is followed by 
taking account of the self-gravity, magnetic field, tur- 
bulence, and gas pressure. The ideal magnetohydro- 
dynamics (MHD) is assumed here for simplicity. The 
barotropic equation of state (EOS) is assumed where 
the gas temperature is 10 K below the critical density 
Per = 2 X 10"^^ g cm-^ (ncr = 5.24 x 10^° cm"^), and 
it increases with the adiabatic index 7 = 7/5 above 
per- This change in temperature reproduces the for- 
mation of the adiabatic core, which corresponds to the 
first core of iLarsonl ([l969). The value of the criti- 
cal density Pcr is taken from the num erical results of 
iMasunaga. Mivama. fc Inutsukal (|1998l ). who studied the 
spherical collapse of molecular cloud cores with radia- 
tion hydrodynamics. Recent multidimensional simula- 
tions using radiation MHD produced results for proto- 
stellar collapse only quantitatively di fferent from simu- 
lations assuming a barotropic EOS (|Commercon et all 
120101 : iTomida et al.Ji2010|). The significant differences 
are restricted within the inside and proximity of the first 
core. 

3. NUMERICAL METHODS 

We calculated the evolution of the cloud cores using the 
AMR code, SFUMATO (Matsumoto 2007). It adopts 
a block-structured grid as the grid of the AMR hierar- 
chy. The total variation diminishing (TVD) cell-centered 
scheme is adopted as the MHD sol ver, with the hyper - 
bolic divergence cleaning method of Dedner et al.l ()2002f ). 
The MHD solver achieves second-order accuracy in space 
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and time. The self-gravity is solved by the multigrid 
method, exhibiting spatial second-order accuracy. The 
numerical fluxes are conserved by using a refluxing proce- 
dure in both the MHD and self-gravity solvers. Periodic 
boundary conditions are imposed. 

The computational domain of [—2Rc^2Rc\^ is initially 
resolved by a uniform grid having 256^ cells. The initial 
resolutions are Ax = 2.8 x 10~^pc, 3.7 x 10~^pc, and 
5.3 X 10"^ pc for / = 1.68, 3.0, and 6.0, respectively. The 
Jeans condition is employed as a refinement criterion; 
blocks are refined when the Jeans length is shorter than 
8 times the ceh width, i.e., (7rc^/Gp)^/^ < 8 Ax (c.f., 
iTruelove et al.lfTQQTl ). where Cp denotes the sound speed, 
and it is a function of density in the barotropic EOS. The 
finest resolution is 7.0 x 10~^ AU for the typical case. 

We followed the collapse beyond the stage in which the 
maximum density exceeds Pmax = 10^ ^po = 10~^gcm~^ 
(^max = 2.62 X 10^^ cm~^) for all the models except for 
the model with fast turbulence and a strong field where 
(a,A^,/) = (0.5,3.0,1.68). In this model, we confirmed 
that the cloud does not undergo collapse even by t = 
15.5tff (= 3.28 X 10^ yr). 

4. RESULTS 
4.1. Less massive cloud cores 
4.1.1. Overview 

Less massive cloud cores, corresponding to models 
with / = 1.68, are examined in this section. All 
the models with / = 1.68 except for the model with 
= (0.5,3.0,1.68) undergo collapse. After the 
collapse begins, the maximum density of the cloud core 
increases rapidly. After it exceeds the critical density 
of the EOS (pcr), the cloud core forms an adiabatic 
core supported against gravity mainly by the thermal 
pressure. The density of the adiabatic core is typically 
p > 10~^^gcm~^, which is approximately 10^ times 
larger than the critical density of the EOS. Th e adia- 
batic core corresponds to the first core of Larson (1969), 
and hereafter we refer to it simply as the first core. Mod- 
els with larger a and/or A4 have longer latency before 
the initiation of the collapse and hence first core forma- 
tion. The first core forms at t = 2.39tff for a model with 
{a,M, f) = (0.1, 1.0, 1.68), and at t = 11.5tff for a model 
with {a,MJ) = (0.25,3.0,1.68). 

Figure [1] shows the evolution of the velocity dispersions 
in the isothermal collapse phase (pmax < Per) for models 
with moderate magnetic fields {a = 0.25). The velocity 
dispersions within the dense reg ion o f p > O.lpmax are 
calculated according to equation (|C10[) . The velocity dis- 
persions decrease in the dense region before the collapse 
begins for the turbulent models (solid lines). When the 
collapse sets in, the velocity dispersion is subsonic even 
for the strong turbulent model {M. = 3.0) as denoted by 
the red solid line. In other words, decay of the turbu- 
lence promotes the collapse in the dense region. This is 
consistent with molecular line observations of the dense 
cores where narrow line widths are obtained. The ve- 
locity dispersion is smaller in the dense region than in 
the whole computational domain (dotted lines) when col- 
lapse sets in. This indicates that the turbulence decays 
in the collapsing dense region selectively, and the other 
region remains turbulent even after the collapse sets in. 

As the collapse proceeds, the velocity dispersion in- 
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Figure 1. Velocity dispersion, (Av), as a function of time in 
the isothermal collapse phase (pmax < Per)- Red, blue, and green 
lines correspond to models with (a, M, f) = (0.25, 3.0, 1.68), (0.25, 
1.0, 1.68), and (0.25, 0.0, 1.68), respectively. Solid hues denote 
velocity dispersion in the dense region of p > O.lpmax, while the 
dotted lines denote velocity dispersion in the whole computational 
domain. Dashed lines is same as the solid lines but for the radial 
velocity, (Avr). Filled circles associated with the solid lines denote 
the stages of pmax = 10^ po (n = 1, 2, • • • , 6). 
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Figure 2. Three-dimensional structures of density and magnetic 
fields for models with / = 1.68. The stage pmax = 10~^gcm~^ 
(rimax = 2.62 X 10-*^^ cm~^) is shown for the collapse models, while 
the final stage of t = 15.5tff is shown for the non-collapse model 
of {a,A4,f) = (0.5,3.0,1.68). Isosurfaces denote the iso-density 
surfaces of p = 3.16 x IQ-^Ogcm-^ (n = 8.28 x lO^cm-^). Tubes 
denote magnetic field lines. The boxes enclose the entire compu- 
tational domain (0.712 pc)^. 

creases and it exceeds the sound speed in the dense re- 
gion. The increase in the velocity dispersion is attributed 
to t he in fall motion as denoted by the dashed lines (see 
eq. |C11| ). The radial infall dominates over the velocity 
dispersion in the stages of Pmax ^ 10~^^gcm~^. Note 
that the initial stage exhibits (Av) /cg 2 for the model 
with A4 = 3.0, and (Av) /cg — 0.7 for the model with 
M = 1.0. These reductions of {Av) are caused by a 
density weighted average in the calculation of {Av). 
Figure [2] shows the cloud structures on the scale of 
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Figure 3. Column density distributions of the entire compu- 
tational domain (0.712 pc)^ at pmax = 10~^gcm~^ (rimax = 
2.62 X 10^^ cm-3) for the collapse models with / = 1.68. For the 
non-collapse model of (o;, A^, /) = (0.5, 3.0, 1.68), the final stage of 
t = 3.28 X 10^ yr is shown. The column densities are shown on a 
logarithmic scale in the x — z planes. 

the whole computational domain, (0.712 pc)^, at the 
stage with Pmax = 10~^gcm~^ for the collapse mod- 
els, and at t = 15.5tff for the non-collapse model with 
(a, A1, /) = (0.5, 3.0, 1.68). The iso-density surfaces in- 
dicate the structures of the cloud cores at the low density 
of p = 3.16 X 10"^°gcm-^ (n = 8.28 x 10^ cm"^), cor- 
responding to the boundary between the cloud core and 
the parent cloud. The models with a high Mach number 
exhibit a density structure highly disturbed by the turbu- 
lence (Fig. [2]a-[2]c), where the initial configuration of the 
cloud core disappears. In models with a moderate Mach 
number, the cloud core is strongly disturbed (Fig.[2]GH2]f). 
The models shown in Figure [2]e and [2]f produce cores that 
are flattened in the plane perpendicular to the magnetic 
field. The models without turbulence produce axisym- 
metric oblate cores flattened in the plane perpendicular 
to the magnetic field (Fig. [2|g'-[2pi). The magnetic field 
lines are highly disturbed by the turbulence in the weak- 
field models (Fig.[2]a and[2]G?). In contrast, the strong- field 
models exhibit almost straight field lines (Fig. [2]c, [2]f, and 

Figure [3] shows the column density distributions of the 
whole computational domain for the models shown in 
Figure [21 clarifying the fine density distribution. The 
cloud cores are shown on a density scale ofn 10^ cm~^, 
corresponding to the C^^O cores. In the model with a 
high Mach number and weak field (Fig.[3]a), the low den- 
sity region is highly disturbed by the turbulence, and 
exhibits a complex structure. The model with a high 
Mach number and strong magnetic field (Fig. [3]c) forms 
a molecular gas sheet, which lies at the top boundary 
of the computational domain. The location of the sheet 
depends on the seed of the random initial velocity field. 
This model shows a filamentary structure aligned parallel 
to the ma gnetic field. A similar structure was reported by 
iPrice fc Bate. (2008>) . who performed an SPH simulation 
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Figure 4. Same as Figure E] but for a (4000 AU)^ region. 
Isosurfaces denote the iso-density surfaces of p = 10"-'^'^ g cm~^ 
(n = 2.62 X 10^ cm-3). The model with (a, M, f) = (0.5, 3.0, 1.68) 
is not shown because it does not undergo collapse and the maxi- 
mum density is lower than the level of the iso-surface. 

taking account of the magnetic field. iNakamura fc Lil 
(|2008l ) also reported such a filamentary structure by per- 
forming grid based MHD simulations including ambipo- 
lar diffusion. Oblate cloud cores form in Figures [3]e, [3]f, 
[3]g', and[3Fi. The cloud cores are more flattened when the 
magnetic field is strong. In the turbulent models, the 
edges of the oblate cloud cores are warped. 

At the stage shown in Figure [3l all the turbulent cloud 
cores move at subsonic speeds. When we define the 
cloud core as the gas denser than the initial ambient gas 
{p > Po/14), the velocities of the baricenter of turbu- 
lent cloud cores range from 0.01 to 0.04 km s~^. Among 
ah the models, the model of (a, M, f) = (0.25, 3.0, 1.68) 
produces the first core at the most distant point (0.39 pc) 
from the initial peak of the cloud core. 

Figure [H shows the density and magnetic field in the 
dense regions of (4000 AU)^ for the collapse models, 
showing iso-density surfaces of p = 10~^^gcm~^ (n = 
2.62 X lO^cm"^). The cloud cores shown in this figure 
correspond to the dense cores observed by dust contin- 
uum emissions. Even on this scale, the weak-field models 
(Figs. [4la and [He) exhibit a disturbed density structure, 
and the magnetic field lines are not aligned. The re- 
maining models exhibit oblate shapes with a minor axis 
parallel to the local mean magnetic field. The magnetic 
field lines are in the configuration of an hourglass. The 
direction of the mean magnetic field depends on the Mach 
number and the magnetic field strength. The direction 
of the local magnetic field in models with a higher Mach 
number and weaker magnetic field is inclined more from 
that of the initial (global) magnetic field. The models of 
Figures [He, [3]f, and [Hg exhibit mean magnetic fields par- 
allel to the 2:-direction, while the models of Figures [H?), 
[4]c, and[4]ti exhibit mean magnetic fields oriented in other 
directions. 

Figure [5] shows the density and the magnetic field on a 
400 AU scale, showing iso-density surfaces of p = 3.16 x 
10"^^gcm-^ (n = 8.28 x 10^ cm"^). On this scale, ah 
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Figure 5. Same as Figure [2]but for a (400 AU)^ region. Isosur- 
faces denote the iso-density surfaces of p = 3.16 x 10"-*^^ g cm~^ 
{n = 8.28 X lO^cm-S). 
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Figure 6. Densities, magnetic fields, and outflows at the out- 
flow formation stages for the model with / = 1.68 in the dense 
region with (40 AU)^. The sta ges are pmax — {cl) 1.08 x 
lO-Sgcm-3, (b) 1.01 X lO-Sgcm-3, (c) 7.55 x lO-^gcm-^, (d) 
1.27X 10~^ gcm-3, (e) 2.16 x 10~^ g cm-^, (/) 1.43 x 10~^ g cm-^, 
(g) 1.04 X 10~^gcm~^. The green isosurface represents a density 
of p = 1.00 X 10-12 gcm-3 (n = 2.62 x lO^^ cm-^). The blue 
isosurface is for the radial velocity of Vr = 0.57 km s"-*^ (vr = 3cs). 
The tubes indicate the magnetic field lines. 

the models except for (a, = (0.1,3.0, 1.68) exhibit 

flat disks in the plane perpendicular to the hourglass- 
shaped magnetic fields. The model with (a,A^,/) = 
(0.1,3.0,1.68) produces a highly warped disk, and the 
magnetic field lines are not aligned even on this scale. 

Figure [6] shows the dense region on a 40 AU scale for 
the outfiow formation stage. Four of the five turbu- 
lent collapse models produce outfiows indicated by the 
blue isosurface of the radial velocity of Vr = 0.57 km s~^ 
{vr = 3cs). The radial velocity is measured from the loca- 
tion of maximum density. The outfiow appears ~ 200 yr 



after the maximum density exceeds the critical density 
of the EOS for all the outfiow models. This epoch cor- 
responds to the stages with 

Pmax ^ 10 gem for 

the model of Figure [6]a, and 

Pmax ^10 g cm for 
the remaining outfiow models. All the outfiows are 
ejected in the direction parallel to the local magnetic 
field. The weak-field models (Figs. [6]a and [6]c) result in 
outfiow directions completely different from the initial 
direction of the magneti c field . This is consistent with 
'Matsumot o fc Tomisakal (|2004l ). who examined outfiow 
formation in collapsing clouds, in which the initial mag- 
netic field and rotation axis were not aligned. The out- 
fiows reproduced here are classified into two types: bipo- 
lar and spiral fiows. The outfiow shown in Figure [6]a is 
a bipolar fiow driven by tightly twisted magnetic fields. 
The bipolar fiows in Figure ^ and [6H are accelerated 
by magnetic fields that are less twisted because of their 
higher strength. In these bipolar fiows, the rotation axis 
of the central first core is parallel to the mean direction 
of the local magnetic field. The other weak-field model 
of Figure [6]c shows a spiral fiow, which is qualitatively 
different from the previous three outfiows, showing mag- 
netic field lines wound in a spiral shape, along which the 
gas is accelerated. The rotation axis of the central first 
core is inclined at a large angle of 35 ° to the direction 
of the local mean magnetic field. The model shown in 
Figure [6]e does not produce an outfiow even at the final 
stage with p = 2.16 x 10~^gcm~^. Because of the high 
magnetic field strength in this model, angular momen- 
tum is transferred mainly by magnetic braking instead 
of by outfiow. The relationship between the outfiow, ro- 
tation, and magnetic field are described in detail in §4.31 

Note that we reproduce the very early phase of outfiow 
formation. At the stages shown in Figure [6l the high- 
est outfiow velocity is found for Figure [6]a, exhibiting a 
maximum value of Vr = 20. Tc^ (= 3.9kms~^). The re- 
maining models produce maximum outfiow velocities of 
^ lOcs (= 2kms~^). For all models, the outfiow velocity 
continues to increase until the end of the calculations. 

Figure [3 shows close-up views of the first cores in the 
column density distributions on a 10 AU scale. All the 
collapse models produce spherical first cores with a ra- 
dius of approximately 1 AU. The masses of the first cores 
are ~ 6 x 10~^ Mq at the stage shown in Figure [3 All 
the first cores are surrounded by disk-shaped envelopes. 
The disk shape can be clearly seen in Figures [7H, [3e, 
[Tjf, and [Tig because these disks are being viewed edge- 
on. The outfiows disturb the disk-shaped envelope near 
the first cores (Fig. [7]a and[71c). Note that both the first 
cores and the disk-shaped envelopes rotate slowly. The 
first cores are supported against gravity by their thermal 
pressure. The disk-shaped envelopes, on the other hand, 
are only partially supported by the centrifugal force; con- 
sequently, the infall velocity dominates over the rotation 
velocity. The disk-shaped en velope resembles the magne- 
tized pseudo disk re ported by lShu fc L i' fl99T). Recently, 
iMellon fc Lil (j2008L l2009f ) also reported that a centrifu- 
gally supported disk cannot be formed in the presence of 
a magnetic field due to magnetic braking. 

4.1.2. Change in shape during collapse 

In order to measure the complexity of the density dis- 
tribution during the collapse, we evaluate the surface-to- 
volume ratio and axis ratios for a given iso-density sur- 
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Figure 7. Column density distributions of central (10 AU)^ re- 
gion for pmax = 10~^ gcm~^. for the models with / = 1.68. The 
column densities are shown on a logarithmic scale in the x — z 
planes. 
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Figure 8. (a) Normalized surface-to- volume ratios and (b) 
axis ratios as functions of density thresholds for the model with 
{a,MJ) = (0.25,1.0,1.68). The different colored lines corre- 
spond to the stages of pmax/po = 10, 10^, lO'^, lO^^, and lO^^ 
(pmax = 10-18, 10-1^, 10-12, 10-9, lO-Sgcm-3). From left to 
right, the solid vertical lines represent the initial density of the 
ambient gas (po/14.0), the initial central density, and the critical 
density of the EOS (pcr)- The dashed vertical lines represent the 
densities at which grid refinement was carried out according to the 
Jeans condition. 

face according to Appendix [B1 The surface-to- volume 
ratio is normalized so that it has a value of unity for 
a sphere. Figure [8] shows surface-to-volume ratios for 
the model with {a,MJ) = (0.25,1.0,1.68) as a fidu- 
cial model. In the early stage with Pmax = 10~^^gcm~^ 
(blue curves), the iso-surfaces of p :^ 10~^^gcm~^ have 
high surface- to- volume ratios, indicating that the cloud 
core has a complex density distribution in the low density 
region. The surface-to- volume ratio vanishes below the 
minimum density. This is because the isosurface vanishes 
there while the volume coincide with that of the whole 
computation box because of the periodic boundary con- 
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Figure 9. (a) Normalized surface-to- volume ratios and (b) axis 
ratios as functions of density thresholds for the models with / = 
1.68 at the stage with pmax = 10-^ gcm-^. The different colored 
lines correspond to the different models. From left to right, the 
solid vertical lines represent the initial density of the ambient gas 
(po/14.0), the initial central density, and the critical density of 
the EOS (pcr). The dashed vertical lines represent the densities 
at which grid refinement was carried out according to the Jeans 
condition. 

dition. In contrast, the iso-surface of p ^ 10~^^gcm~^ 
exhibits a low surface- to- volume ratio, and as is consider- 
ably longer than a2 . This indicates that the collapsing re- 
gion has a smooth triaxial shape at this stage. Note that 
the vertical dashed lines denote the densities at which 
the grid refinements are performed. At low densities, the 
surface-to- volume ratio is seen to make a transition from 
a high to a low value. Since this occurs at a density 
less than that of the first refinement, it cannot be at- 
tributed to the grid refinement process. At the stage with 
Pmax = 10~^^gcm~^ (green curves), the two axis ratios 
have an identical value of 3 — 4, indicating that the col- 
lapsing region of the cloud core is an almost axisymmetric 
disk. Owing to its fiat shape, the surface-to-volume ra- 
tio is slightly higher than unity for p ~ 10~^^gcm~^. At 
the stage with Pmax = 10~^^gcm~^ (yellow curves), the 
surface-to-volume ratio reaches 1.8 at p 10~^^gcm~^, 
and the axis ratios reach 7 — 8, indicating that the flat- 
ness of the disk increases during the collapse. At the 
stage with Pmax = 10~^gcm~^ (red curves), a spher- 
ical first core forms, and the surface- to- volume ratio 
and axis ratios become unity at the high density of 
p > 10~^^gcm~^. The surface-to- volume ratio and the 
axis ratios of the disk-shaped envelope reach 2 and 10, 
respectively, at p 10~^^gcm~^ owing to an increase 
in the fiatness. At the outfiow formation stages (purple 
curves), the surface-to- volume ratio takes a high value of 
2.5 at p 10~^^ gcm~^, reflecting the fact that the out- 
flow disturbs the disk-shaped envelope around the flrst 
core. The surface- to- volume ratio and the axis ratios re- 
main unity in the dense region with p > 10~^^gcm~^, 
indicating that the flrst core is spherical even at the out- 
flow formation stage. Moreover, at the low density re- 
gion with p 10~^^gcm~^, the surface-to-volume ratio 
remains high even after flrst core formation, indicating 
that the periphery of the cloud core remains turbulent. 

Figure [9] compares the surface-to- volume ratio and the 
axis ratios for the four models with / = 1.68. For a 
low density of p < 10~^^gcm~^, the turbulent mod- 
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els exhibit high surface-to- volume ratios. The model 
with a high Mach number {A4 = 3) exhibits the high- 
est value (red curve), while the model without turbu- 
lence exhibits a low ratio (yellow curve). Comparing the 
models with = 1, the weak-field model (blue curve) 
has a higher surface-to- volume ratio than the strong-field 
model (green curve), implying that disturbance by tur- 
bulent fiow is considerably suppressed by the magnetic 
field. 

In all the models, the axis ratios increase with density 
in the range of 10~^^gcm~^ ^ P ^ 10~^^gcm~^. In 
this range, the mean axis ratios (a2 + as)/{2ai) tend to 
increase with the initial magnetic filed strength a, and 
decrease with the initial Mach number A4. This implies 
that the magnetic field increases the degree of anisotropy 
and the turbulence increases the effective sound speed. 
This tendency is apparent at a relatively low density of 
g ^^-3 < p < 10-15 g cm-^. 

At p ^ 10~^^gcm~^, the model with {a,M,f) = 
(0.1,1.0,1.68) exhibits a prolate shape (blue curves), 
with the major axis being considerably longer than the 
other axes; the axis ratio as/a2 has a maximum of 6.2 
at p = 3.8 X 10~^^gcm~^. For all the other models, 
the axes as and a2 are comparable {as/a2 < 2) for 
p > 10~^^gcm~^, indicating an oblate shape. For a 
higher density of p > 10~^^gcm~^, all the models ex- 
hibit spherical shapes because of the presence of the first 
cores. 

4.2. Dependence on mass 

We examined two additional models with / = 3.0 and 
6.0 in order to investigate the dependence on cloud mass. 
We refer to these models as massive cloud cores. The 
parameters of the initial magnetic field and the initial 
Mach number of turbulence are a = 0.1 and A4 = 3.0, 
respectively. It takes shorter time to form first core for 
the more massive core; the first core forms at t = 5.24tff 
for a model with (a,A^,/) = (0.1,3.0,1.68), and at t = 
1.29tff for a model with {a,M, f) = (0.1, 3.0, 6.0). 

Figures [10] and [11] compare the three models on three 
spatial scales. The low density regions are highly dis- 
turbed by the turbulence in all the models (left column 
in Figs. [To] and [TT]). On the intermediate scale (the mid- 
dle column), the clouds take the shape of filaments for 
the massive models. The most massive model produces 
a long thin filament as shown in Figures {TOh and {TTh . 
On the small scale, all the models produce spherical first 
cores (right column in Fig. [TT]) . The first cores of the 
massive clouds are embedded in the filament, while that 
of the less massive cloud is surrounded by the disk (right 
column in Figs. [TO]) . 

All the models exhibit outfiows as shown by the blue 
iso- velocity surfaces in Figure [TO] but their shapes are 
different. The outfiow shown in Figure [TOb is classified 
as bipolar fiow, and it is ejected in the direction per- 
pendicular to the plane of the disk. The disk-outflow 
system is roughly axisymmetric although its axis is in- 
clined considerably with respect to the initial direction 
of the magnetic field (the vertical direction). In the mas- 
sive models with / = 3.0 and 6.0 (Figs. ^ and [IS), 
the outflow is classified as spiral flow, and it is ejected 
in the plane perpendicular to the filament, and is not 
collimated. Similar outfiow appears for the model with 
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Figure 10. Three-dimensional structures of density, magnetic 
fields, and outfiows for the models with (a, A4, f) = (0.1, 3.0, 1.68), 
(0.1,3.0,3.0), and (0.1,3.0,6.0) from top to bottom. From left to 
right, the plotted areas are the entire computation box, (400 AU)^, 
and (40 AU)^. From left to right, the isosurfaces denote the iso- 
density surfaces oi p = 3.16 x 10~^° gcm~^, 3.16 x 10~^^ gcm~^, 
and 1.00 x lO"-*^"^ g cm~^. Tubes denote magnetic field lines. The 
blue isosurface is for a radial velocity of Vr = 0.57 km s"-*^ (vr = 
3cs). The stages are pmax = 4.53 x 10~^gcm~^ for panel /, and 
Pmax = 10~^ gcm~^ for the rest of the panels. 
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Figure 11. Column density distributions at the stage with 
Pmax = 10~^gcm~^ for the models with (a, A1,/) = 
(0.1,3.0,1.68), (0.1,3.0,3.0), and (0.1,3.0,6.0) from top to bot- 
tom. The column densities are shown on a logarithmic scale in the 
X — z planes. From left to right, the plotted regions are the whole 
computational domains, (200 AU)^, (10 AU)^. 

{a,MJ) = (0.1,1.0,1.68) as shown in Figure [6]c. The 
configuration of the magnetic field depends on the type of 
outfiow present. In the less massive model with / = 1.68, 
the magnetic field lines are twisted and highly wound 
along the bipolar directions, implying that the outfiow 
is accelerated by the magnetic pressure. In the massive 
models, the magnetic field lines exhibit a spiral configu- 
ration, which is attributed to misalignment between the 
rotation axis and the local mean magnetic field. This 
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Figure 12. (a) Normalized surface-to- volume ratios and (b) axis 
ratios as functions of density thresholds for the models with A4 = 3 
and a = 0.1. From left to right, the solid vertical lines represent 
the initial density of the ambient gas (po/14.0), the initial central 
density, and the critical density of the EOS (pcr)- The dashed 
vertical lines represent the densities at which grid refinement was 
carried out according to the Jeans condition. 
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Figure 13. Distribution of axis ratio for the isothermal range 
(po < P < Per) at the stage with pmax = 10~^gcm~^ for all the 
collapse models with turbulence {JV[ / 0). The open circles denote 
the axis ratios at p = pQ. The filled circles denote the axis ratios 
at p = lO^po = 1,2, •••,6). The dashed lines indicate the 
boundary between prolate, triaxial, and oblate shapes. 

morphology implies that the outflow is accelerated by 
magnetic tension. 

Figure [12] compares the surface-to- volume ratio and 
axis ratios for the three models. Even for the same 
initial Mach number, low-density regions with p ~ 
10~^^gcm~^ are more disturbed in the more mas- 
sive model, showing a higher surface-to-volume ratio 
(Fig. [T2la). The less massive model produces a disk- 
shaped envelope in the range 10~^^gcm~^ ^ P ^ 
10~^^gcm~^, indicated by the identical values of the 
two axis ratios (red curves in Fig. [T2b). In contrast, 
the massive clouds produce a filamentary envelope. The 
axis ratio a^/ai reaches a maximum value of 80 for 
the most massive model (blue curves in Fig. [T2b). For 
p > 10~^^gcm~^, the surface-to- volume ratio and the 
axis ratios are almost unity, irrespective of the cloud 
masses because of the spherical first cores. 

Figure [13] shows the distribution of the axis ratios 
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Figure 14. Energy ratios (a) {E^^^ + £;th)/|^grav | and (6) 
^kin/^th cis functions of the maximum density pmax in the isother- 
mal collapse phase for the models with = 3 and ol = 0.1. Solid, 
dotted, and dashed lines correspond to models with / = 1.68, 3.0, 
and 6.0, respectively. Gray lines in the panel (6) denote the ratios 
of the kinetic energy of the radial velocity to the thermal energy 

^kin,rad/ ^th- 



of the isothermal envelopes {po < p < Pcr)- Based 
on the axis ratios, the shapes of the envelopes are di- 
vided into three cat egories: prolate, triaxial, and oblate 
(|Gammie et al.l [2003 ). In all the collapse models, the 
shape anisotropy increases (ai/aa decreases) with the 
density threshold for p > lOpo (= 10"^^gcm"^). Ah the 
massive cloud cores with / = 3.0 and 6.0 have prolate 
shapes, with a2/as and ai/as decreasing with increas- 
ing density threshold. Four of the five less massive cloud 
cores with / = 1.68 produce oblate cores for p > lOOpo- 
In addition, the evolutionary tracks of the axis ratios for 
the dense regions {p > Pmax/10) exhibit loci similar to 
those shown in Figure [T3l indicating that the shapes of 
the envelopes reflect the history of the collapse. In sum- 
mary, massive cloud cores tend to have a prolate shape, 
and less massive cores an oblate shape. 

In order to investigate the deformation of the 
cloud cores, we consider the kinetic, thermal, and 
gravitational energy of the dense regions during the 
collapse. Figure [Mb shows (£^kin + £^th)/|^grav| 
in the dense region with p > O.lpmax during 
the isothermal collapse phase {po < Pmax ^ 
per)- These energies are estimated as E]^ 

and |i?grav| = (l/87rG)/^>o,^^_Jg, 
in Figure [Mb . the values of the energy ratio (£^kin + 
£^th)/|£^grav| converge to ~ 4 by the early stage with 
Pmax — 10~^^gcm~^, and remain roughly constant 
during the isothermal collapse. This convergence is 
explained by the virial theorem, which predicts that 
(£;kin + £^th + £^mag/2)/|£^grav| - 1/2. The ovcrestimatiou 
in Figure [Mb is attributed to the underestimation of the 
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gravitational energy, since we ignored the contribution of 
the low density regions. 

Figure [T4b shows the energy ratio, £^kin/^th, during 
the isothermal collapse. For Pmax ^ 10~^^gcm~^, the 
kinetic energy exceeds the thermal energy for the most 
massive model with / = 6 (thin dashed curve), and vice 
verse for the less massive model with / = 1.68 (thin solid 
curve). The massive cloud core is supported mainly by 
the turbulence, and it collapses to form a filament. In 
contrast, the oblate shape in the less massive model is 
attributed to thermal pressure support. 

In the early phase, the energy ratio £^kin/^th decreases 
for all the models. However, they begin to increase at 
Pmax = 10~^^gcm~^ for the model with / = 1.68 and 
at 10~^^gcm~^ for the model with / = 3. This in- 
crease is attributed to an increase in the infall energy, 
as indicated by the gray curves, which denote the en- 
ergy ratio ^kin,rad/^th, where ^kin ,rad is the kinetic en- 
ergy associated only with the radial velocity, and is de- 
fined by ^kin,rad = (1/2) S p>Q,lp^^^ P^l^V . Yoi Pmax > 

10~^^ gcm~^, the energy of infall motion dominates over 
the turbulent energy for the models with / = 1.68 and 
3.0. 

4.3. Outflow, rotation, and magnetic field 

As shown in ^4.11 and ^4.2[ the models examined here 
exhibit either bipolar or spiral outflow. Examples of 
bipolar flows are shown in Figures [6]a, [6j), and [^i, and 
spiral flows in Figures [6]c, [TQ]f and [TOt . The bipolar out- 
flows tend to be associated with disk-shaped envelopes, 
while the spiral flows tend to have filamentary envelopes. 

In order to investigate the relationship between the 
rotation and magnetic field, we calculate the angular 
momentum (J) and the mean magnetic field {B) for 
p ^ O.lpmax during the collapse (see Appendix IC]) . The 
angular momentum can be decomposed into two compo- 
nents parallel and perpendicular to the mean magnetic 
field, defined by 

J|I = |J-B|/|B|, (3) 

and 

= |J X B\/\B\. (4) 

It should be pointed out that the region of interest 
p > O.ljOmax changes temporally. Therefore, changes in 
J and B do not indicate temporal changes in the angular 
momentum and magnetic field of the fixed regions, but 
rather they trace the angular momentum and magnetic 
field in the collapsing region. 

Figure [l5] shows J/M^, Jy/M^, and J^/M^ as 
functions of Pmax for the model with (a,A^,/) = 
(0.25,1.0,1.68), the representative model of a bipolar 
outflow. The ordinate J/M^ is nearly equal to the spe- 
cific angular momentum per unit mass j/M, which re- 
mains constant as the disk-like cloud collapses in the 
cylin drical radial direction (see e.g., i Matsumoto et al.l 
|1997[ ). For a slowly rotating non-magnetized cloud, 
J/M^ increases slightly with density because of a 
slight spin-up due to sp herical collapse (see Fig. 4 in 
IMatsumoto fc Tomisakafl2004J ) . 

In isothermal collapse (pmax < Per), J/M'^ decreases 
slightly with a considerable amount of oscillation. The 
slight decrease in the J/M^ is attributed to magnetic 
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Figure 15. Normalized angular momentum ( J/M^ [cs/47rG]) 
within the dense region with p > O.lpmax as a function of 
the maximum density (pmax) for the model with (a,A4,f) = 
(0.25,1.0,1.68). The angular momentum J/M'^ is plotted as a 
thin solid curve, while the parallel and perpendicular components 
of angular momentum with respective to the local magnetic field, 
J||/M^ and J^/M^, are plotted as thin dashed and thin dotted 
curves, respectively. The thick solid curve denotes the angle Oqi 
between the vectors oi the angular momentum J and the mean 
local magnetic field B for the region with p > O.lpmax, while the 
thick dashed curve denotes the angle 6cr, which is same as ^oi but 
for the region of p > pcr- Both angles are restricted to within [0, 90] 
degrees. The vertical line denotes the critical density of the EOS, 
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Figure 16. Loci of the directions of the mean magnetic field 
±S (red line), the angular momentum ±J (green line), and the 
minor axis of the density structure (blue line) in the region with 
p > O.lpmax for the model with (a, M, f) = (0.25, 1.0, 1.68). Digits 
n with filled circles denote the stages of pmax = 10^ po- Black filled 
circles denote the final stage. The three large circles denote the 
relationships {$1 + 6'2)i/2 = 30°, 60°, and 90°. 

braking. In the early phase with Pmax ^ 10~^^gcm~^, 
the perpendicular component, Jx, is larger than the par- 
allel component, Jy, and exhibits a large value of ^01, 
which is the angle between J and 5, and it is restricted 
to values from O to 90°. The angle ^01 therefore becomes 
when J and B are either parallel or anti-parallel. In 
the range of 10~^^gcm~^ < pmax ^ 10~^^gcm~^, J± is 
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Figure 17. Same as figure [TSl but for the model with (a, /) = 
(0.1,3.0,1.68). 

smaller than Jy , and ^oi also becomes smaller. When the 
EOS becomes adiabatic (pmax = Per), is larger than 
J||. Also plotted in Figure [15] is another angle ^cr, which 
is the angle between the angular momentum and the 
magnetic filed for the region with p > per- After the EOS 
becomes adiabatic (pmax ^ Per), the angle Ocr remains 
constant. In contrast, ^oi decreases up to ~ 5° by the fi- 
nal stage, indicating that the rotation axis is aligned with 
the magnetic field in the dense region. Moreover, the an- 
gular momentum J/^M decreases. The perpendicular 
component decreases selectively by more than an order 
of magnitude. This indicates that the magnetic braking 
is more effective against the perpendicular component of 
angular momentum with respective to the local magnetic 
field than the parallel component . Such selective mag- 
netic braking is also reported in lMatsumoto fc Tomisakal 
([2004). 

Figure [16] shows the loci of the directions of the mean 
magnetic field the angular momentum J, and the 
minor axis of the density structure (the normal vector of 
the disk) in the region with p > O.lpmax for the model 
with {a,MJ) = (0.25, 1.0, 1.68). Each vector is plotted 
in the two-dimensional plane as 



8iictdin{Vxy/Vz) ( 



xy 



(5) 



where V = {Vx^Vy^Vz) represents a given vector and 
Vxy = {Vx + )^^^- The distance from the origin rep- 
resents the angle between a given vector and the z-axis. 
For example, a vector parallel to the x-axis is plotted 
at (90°, 0). The direction of the magnetic field and the 
minor axis of the density structure are aligned during 
the collapse, suggesting that the cloud collapses along 
the magnetic field lines, and the disk is oriented always 
perpendicular to the local magnetic field. The direction 
of the angular momentum drifts considerably over the 
— plane, while it begins to be aligned with the 
magnetic field and the minor axis at Pmax/po ~ 10^ 
(Pmax = 10-11 gcm-3). Finally, for Pm.^/po > 10^^ 
(Pmax = lO-^gcm"^), all three vectors are aligned at 
the point {Ox.Oy) (-40°, -20°). This alignment pro- 
duces the bipolar outflow. 
Figure [17] is same as Figure [15] but for the model with 




Figure 18. Same as figure fT6] but for the model with (a, A^, /) = 
(0.1,3.0,1.68). The loci are shown for pmax / po > 10, because they 
are noisy below that stage. 

with (a,A^,/) = (0.1,3.0,1.68), where the most promi- 
nent bipolar outflow appears. In the early stage with 
Pmax ^ 10~^^gcm~^, the angular momentum strongly 
oscillates due to the strong turbulence and the weak 
magnetic field. After the collapse proceeds, J/M^ re- 
mains roughly constant in the isothermal collapse phase 
(Pmax ^ per)- In this phasc, the direction of the an- 
gular momentum changes considerably as shown in Fig- 
ure [18] The weak magnetic field changes its direction sig- 
nificantly due to precession. The disk normal follows the 
magnetic field, so that the plane of the disk remains per- 
pendicular to the magnetic field. In the adiabatic phase 
(Pmax > Per), J/M'^ dcccascs with increasing density 
as shown in Figure [T7] The perpendicular component 
J±/M'^ decreases drastically as in the previous model. 
This model therefore results in alignment of the rotation 
axis with the local magnetic field, with Oqi in the 
final stage. This alignment is also shown in Figure [181 
where the directions of the magnetic field, the angular 
momentum, and the disk normal converge at the point 
i^x^^y) — (80°, 30°), that is almost perpendicular to the 
z axis. The angular momentum J /M'^ is about five times 
larger than that for the previous model throughout its 
evolution (compare FigsITS] and [TTj) . and this model ex- 
hibits strong outflow. The large angular momentum is 
attributed to both the initial strong turbulence and the 
weak magnetic field. 

Figure [19] shows the evolution of angular momentum 
for the model with (a,A1,/) = (0.1,3.0,3.0), a typical 
model showing a spiral flow. In this model, the angular 
momentum and magnetic field are not aligned through- 
out the evolution. The angle ^oi increases and exhibits 
large oscillations. The angular momentum J/M^ re- 
mains roughly constant in the isothermal collapse phase 
(Pmax < Per), whilc it dccrcascs in the adiabatic phase 
(Pmax > Per)- The pcrpeudicular component J±_ is con- 
siderably larger than the parallel component Jy in the 
adiabatic phase, producing large ^oi and Ocv values. 

The magnetic field, the angular momentum, and the 
minor axis evolve completely differently in the spiral 
models as shown in Figure [20] The locus of J moves in 
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Figure 19. Same as figure [TSl but for the model with (a, /) 
(0.1,3.0,3.0). 




Figure 20. Same as figure [TgI but for the model with (a,A4, f) = 
(0.1,3.0,3.0). 

the lower right direction with a large precession. The lo- 
cus of B moves diagonally from lower left to upper right, 
indicating that the orientation of B rotates around J. 
The point B traverses the 6x — Oy plane three times, in- 
dicating that it makes 1.5 revolutions. The minor axis 
follows the direction of the magnetic field in the stage 
with Pmax < lO^Vo (pmax < IQ-^gcm"^), while it fol- 
lows the direction of the angular momentum in the stage 
with Pmax > lO^Vo (Pmax > 10"^gcm"^). In the final 
stage, the angular momentum and the minor axis are 
aligned at the point {Oa,,Oy) (-30°, 80°). This indi- 
cates that the filaments shown in Figure [TOl e and Fig- 
ure [TTI e are perpendicular to the magnetic fields, while 
the first core shown in Figure [TTIe is slightly oblate, with 
its minor axis oriented parallel to the rotation axis. 

5. DISCUSSION 

5.1. Shape of cloud core 

Cloud core shape is investigated quantitatively by 
means of the surface-to-volume ratio and axis ratios 
(Figs, in [HI and [T2]). In all the models with turbulence. 



a cloud core has a complex density distribution before 
it begins to collapse. In contrast, the density distri- 
bution is relatively smooth in a collapsing cloud core, 
which is partly due to decay of turbulence before col- 
lapse. Strong turbulence prohibits cores from collapsing. 
In other words, the decay of turbulence promotes collapse 
in particular in a less massive core (Fig.[T]). A dense core 
begins to collapse when the turbulence decays to a cer- 
tain level depending on the core mass. A less massive 
core does not collapse until the turbulence decays to a 
very low level. Accordingly the collapsing cores tend to 
have a smooth density distribution. We also note that 
it takes a longer time before collapse if the initial tur- 
bulence is stronger. This is one of the reasons why the 
density distribution is relatively smooth in the first cores 
formed in our simulations. 

The above argument, however, does not hold for a mas- 
sive cloud core, which collapses rapidly before decay of 
turbulence (Fig. [T4|) . Even when a core forms before 
the decay of turbulence, it has a relatively smooth den- 
sity distribution. We suppose that this is due to shrink 
of the core. The diameter of a collapsing core is about 
the Jeans length, which decreases in proportion to p"^/^ 
during the runaway collapse if the gas is isothermal. Tur- 
bulence supports a core against collapse only when the 
typical wavelength is shorter than the diameter of the 
cloud. Turbulence has a smaller power at a shorter wave- 
length in our initial model as well as in the Kolmogorov 
spectrum. Accordingly the effective level of turbulence 
decreases along with the shrink of the core, unless the 
wavelength of the turbulence shortens in proportion to 
the Jeans length. The wavelength of turbulence may be 
shortened a little by compression and advection, but not 
much. Remember that the core mass decreases during 
the runaway collapse. This means that the advection is 
slower than the collapse. The above mentioned idea is 
supported by the fact that the turbulence is still strong in 
the low density region even after the collapse (Fig. [12]). 

Consequently, the turbulence introduces only smooth 
motion, e.g., rotation and shear, into the collapsing re- 
gion. Our results are in agreement with observations of 
starless dense cores with molecular line emissions (e.g., 
iGoodman eFaIlll998l : iTafalla et ani2002[ ). which suggest 
that turbulent motion decreases toward the center of the 
cloud core and coherent motion remains at the central 
part. For star- forming cores, the internal structure may 
be produced by the activation of young stars, e.g., out- 
fiows from young stars and the orb ital motion of a mul- 
tiple system fc.f.. IWang et al.ll2010[ ). 

The shape of the collapsing region depends mainly on 
the mass of the cloud core, and the orientation is con- 
trolled by the magnetic fields. In the less massive cloud 
core, the collapsing region assumes an oblate shape, 
whereas in the massive cloud cores, it assumes a pro- 
late shape. In all cases, the minor axis is parallel to the 
local magnetic field and the shape anisotropy increases 
with density in the isothermal envelope. 

The deformation during the collapse can be explained 
by the energy ratio as shown in Figure [TH This analysis 
indicates that the pressure exerts an isotropic effect, and 
it impedes deformation of filaments in the less massive 
cloud cores. In other words, the massive cloud core is 
subject to few isotropic effects, and it becomes deformed 
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into a filament. Thi s is consis t ent w ith the results of a 
classical analysis by iLin et aP ()1965l ) , who investigated 
deformation during dust cloud collapse. 

There have been many simulation-based studies on 
large-scale turbulence, in which the shap es of clum ps and 
cloud cores are have been considered. iLi et al.l (|2004) 
performed 512^ simulations with magnetic fields, and re- 
ported tha t the majority of cores are prola te or triaxial 
i n sha pe. lOffner Krumholzl (j2QQ9[ ) and lOffner et al] 
poos) performed large-scale AMR simulations without 
magnetic fields, and showed that the cores are predom- 
inantly triaxial. Our simulations demonstrate that the 
shape of the cloud core evolves with increasing density 
(Fig. [13]), though a massive core tends to be pr olate and a 
less massive core tends to be oblate. lOffner fc Krumholzl 
(I2OO9) also suggested that the final shape is somewhat 
sensitive to a cutoff density. 

We found that the minor axis tends to be aligned with 
the local magnetic field irrespective of the shape of the 
cloud cores (oblate or prolate) in the isothermal collapse 
phase. A similar tendency was also found by Li et al. 
([2004) , who reported a weak correlation between the mi- 
nor axis of the cloud core and the local magnetic field. 
Our simulations indicate that the minor axis is aligned 
rapidly with the local magnetic field during the isother- 
mal collapse phase (e.g.. Fig. [T6|). Moreover, in the weak 
magnetic field models, the minor axis changes orienta- 
tion during the collapse, following the direction of the 
magnetic field (see Figs. [181 and [20]). Therefore, the mi- 
nor axis is always oriented parallel to the local magnetic 
field. 

Hourglass-shaped magnetic fields are reproduced in the 
present simulations, as shown in Figures [4] and [3 which 
are on 4000 AU and 400 AU scales, respectively. The 
hourglass shapes are prominent in the models with mag- 
netic fields stronger than or equal to the fiducial strength 
{B > 20 /iG at n ~ lO^cm"^). Observations of po- 
larization have revealed such hourglass-sh aped magnetic 
fields in both high-mass (Schl euninglll998[ ) and low-mass 
(|Sugitani et a l. 2010) star-forming cores. The spatial 
scales involved are ~ 1 pc, which are larger t han the scale 
of the present study. A notable example is iGirart et al.l 
()2006[ ). who resolved the magnetic field in the low-mass 
core, NGC 1333 IRAS 4A, on a scale of a few hundred 
AU. They revealed hourglass-shaped magnetic fields per- 
pendicular to the elongated envelope. However, NGC 
1333 IRAS 4A is a binary protostar system. The possi- 
bility of binary formation is discussed in 



(a) Bipolar flow 



(b) Spiral flow 



5.2. Outflows 

The models here exhibit either bipolar or spiral 
outfiows, as illustrated in Figure [2TJ The former 
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120091 ). The results of our simulations demonstrate that 
bipolar fiows are also produced in turbulent cloud cores. 

The bipolar flows are divided into two subtypes ac- 
cording to the configuration of the associated magnetic 
field. The outfiows in the moderate field models have 
a poloidal magnetic field larger than the toroidal field, 
as shown in Figures [3) and [6]ri. Another subtype is re- 
ported in Figure [6]a, where the toroidal magnetic field 




Figure 21. Schematic diagram of two types of outflows: (a) 
bipolar flow and (b) spiral flow. The surfaces represent iso-density 
surfaces, the tubes denote the magnetic fleld lines. The arrows 
indicate the direction of the outflows. 

dominates over the poloidal field. The former out- 
fiow is driven by the magneto- centrifugal mechanism 
(c.f. , 'Blandford & Payne"1982'; 'Pudr itz fc Normanl[T986l ) , 
and the latter by the magnetic pressure gradient of the 
toroidal field co mponent, categorized as I- type fiow in 
iTomisakal (|2002 f). 

The spiral fiow is a previously unreported type of 
outfiow. The magnetic field lines are wound up by 
the rotation of the axes oriented perpendicular to the 
fiel d. A similar ni agnetic field morphology was reported 
by iMachida et al.l f|2006) for inclined rotators (see their 
Fig. 13). However, they did not confirm this type of out- 
fiow in their simulations because they could not follow 
the evolution for a sufliciently long period. 

It is not yet known whether such spiral flows are stable 
over periods of more than several thousand years. When 
the angular momentum perpendicular to the magnetic 
field is released by the outfiow, the outfiow may trans- 
form from spiral to bipolar. In order to confirm this pos- 
sibility, we need to employ the sink particle method to 
simulate longer timescales at a reasonable computational 
cost. 

5.3. Probability of fragmentation 

Figure [22] displays the evolution of magnetic fiux- 
spin relations for the models with M. > 0. The mag- 
netic fi ux-spin relations were proposed bv IMachida et all 
as a means of investigating fragmentation 
during cloud collapse. The abscissa and ordinate de- 
note non-dimensional parameters of magnetic fiux b = 



|S|/(87rc^p)V2 and spin 



|l^|/(47rGp)i/2, respec- 



tively, for the dense region with p > O.lpmax (see Ap- 
pendix [Cj. This diagram predicts fragmentation driven 
by rotation. 

The solid curves in Figure [221 are loci in the isothermal 
phase (pmax < per) for all the turbulent models. Starting 
from the initial condition (open circles), all the loci drift 
in the b — uo plane, until they reach the convergence curve 
denoted by the black curve. This convergence was also 
seen for the cas e of collapse in the absence of turbulence 
f Machida et a l. 2005a. b). They summarized the condi- 
tions for fragmentation as follows. If the locus reaches 
the horizontal part of the convergence curve (rotation 
dominant) during the isothermal collapse, the cloud core 
fragments. If it reaches the vertical part (magnetic field 
dominant), it does not. All the models examined here 
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Figure 22. Loci in the magnetic flux-spin relation plane for 
the models with > 0. The abscissa and ordinate denote non- 
dimensional parameters of magnetic flux (6) and spin (uu), respec- 
tively. Open circles, fllled circles, and diamonds indicate the stages 
oi p = po (initial stage), pcr, and p = 10 ■'^■'^po- Loci with con- 
tinuous curves are for the isothermal phase (pmax < Per), and 
loci with dotted curves are for the adiabatic phase (pmax > Per)- 
The quarter oval curve indicates the convergence curve, given by 
(6/0.36)2 ^ (^/o.2)2 = 1. 

reach the vertical part of the convergence curve in the 
isothermal collapse phase, and they do not undergo frag- 
mentation. This indicates that turbulent cloud cores 
have insufficient angular momentum to fragment. 

Another possibility for fragmentation is turbulent frag- 
mentation. In the massive models, the turbulence ex- 
ceeds the thermal pressure as a supporting force against 
gravity. These models produce thin filaments. Each fil- 
ament produces only one first core in the present simu- 
lations, because of the very small time steps in the adia- 
batic phase. However, the formed filaments are very thin 
and tend to undergo additional fragmentation in the later 
stages. Such fragmentation could be reproduced if sink 
particles were introduced into the simulations. 

After the maximum density exceeds the critical den- 
sity Per, the loci traverse the convergence curve and they 
move in the leftward direction (dotted curves). This left- 
ward movement is attributed to an increase in the sound 
speed Cp in the adiabatic phase (pmax > Per)- Some 
loci move close to the horizontal part of the convergence 
curve, where rotation support exceeds magnetic field sup- 
port. Based on these models, a protoplanetary disk may 
fragment in the later stages if it is sufficiently massive 
compared to its central star. 

6. SUMMARY 

The collapse of turbulent magnetized cloud cores is in- 
vestigated by AMR simulations, resolving both the cloud 
core and the first core. 

The cloud core has a complex density distribution in 
the low density region corresponding to its boundary 
with the parent cloud, because of disturbances due to 
turbulence. After the collapse begins, the density dis- 
tribution becomes smooth in the collapsing region. This 
indicates that the collapse dilutes the tiny fluctuations 



caused by the turbulence. Even after formation of the 
first core, the edge of the cloud core remains turbulent. 

The shape anisotropy of the collapsing region increases 
during isothermal collapse and depends mainly on the 
mass. When a cloud core is less massive (/ = 1.68), the 
dense region of the cloud core is oblate, with its minor 
axis parallel to the local magnetic field. The collapsing 
region is threaded by an hourglass-shaped distribution 
of magnetic field lines on a scale of < 1000 AU. When a 
cloud core is massive (/ = 3.0 and 6.0), the dense region 
of the cloud core is prolate, with the minor axis parallel 
to the local magnetic field. The extremely massive cloud 
(/ = 6.0) produces a very thin filament, which may later 
fragment. 

In all cases, the orientation of the cloud core is con- 
trolled by the magnetic field. The minor axis of the 
collapsing region changes direction during the collapse, 
following the direction of the local magnetic field irre- 
spective of the cloud core shape. Each model produces 
a spherical first core, in which the density is higher than 
10~^^gcm~^. The first core is embedded in the in- 
falling envelope. 

The dependence of the shape on the mass can ex- 
plained by the energy ratio between the thermal pres- 
sure and the turbulence. When the energy due to ther- 
mal pressure exceeds that due to turbulence, the cloud 
core becomes oblate in the early stages of collapse. If 
the opposite is the case, then the cloud core assumes a 
filamentary shape. 

We found two types of outflows: bipolar and spiral 
flows. Bipolar flow is associated with the disk-shaped 
envelope, while the spiral flow is associated with the fil- 
amentary envelope. Bipolar flow tends to occur in less 
massive cloud cores, and the rotation axis, magnetic field, 
and the disk normal become aligned in the proximity of 
the first core. The disk-outflow system is inclined com- 
pletely from the global magnetic field for the weak mag- 
netized cloud core {B = 7/iG). Even for the moderate 
field models {B = 19/iG), the outflow direction diverges 
considerably from that of the global magnetic field by 
an angle of 40 — 50°. For the strong field model 
{B = 37/iG), the disk envelope is aligned with its mi- 
nor axis along the direction of the global magnetic field, 
while the outflow is not produced. The angular momen- 
tum perpendicular to the local magnetic field is selec- 
tively reduced in the dense region by magnetic braking 
and the outflow during the adiabatic phase in the bipolar 
flow models. The spiral flow tends to be associated with 
massive cloud cores, and the rotation axis is not aligned 
with the magnetic field in the envelope surrounding the 
first core. 

Numerical computations were carried out on a Cray 
XT4 at the Center for Computational Astrophysics 
(CfCA) at the National Astronomical Observatory of 
Japan. This research was supported in part by a 
Grant-in-Aid for Scientific Research (C) 20540238 and 
(B) 22340040 from the Ministry of Education, Culture, 
Sports, Science and Technology, Japan. 
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APPENDIX 

A. ENERGIES OF THE CRITICAL BONNER-EBERT SPHERE 

Th e gravitational potential ^ of the Bonner-Ebert sphere is obtained by solving the equations (see iChandrasekharl 
[1931) 

imposing the boundary conditions of k;(0) =0 and ?/(0) = 0, where the non-dimensional variables ^ and w are defined 

by 

? = ^, (A3) 
w=% = -hi(^Y (A4) 

ci \poJ 

The critical Bonner-Ebert sphere has a maximum radius of ^max = 6.45 (rmax = 6.45a). The central gravitational 
potential is equal to zero, <l>(0) = 0, owing to the boundary conditions. Considering an isolated critical Bonner-Ebert 
sphere, the gravitational potential in r < rmax is given by 

*(r) = $(r)-$(r„,ax)-G^^^^^, (A6) 

'^max 

where M(^) denotes the mass inside the radius ^ = r/a, 

M(0 = 47r / pr^dv = 47ra^po^(0^^- (A7) 
Jo 

The gravitational potential is therefore obtained as 

^(r) = Cl [^(0 - ^(^max) - ^maxy(^max)] • (A8) 

The gravitational energy is given by 

^grav = J / P'^dV (A9) 

= -352pocla^ (All) 

The thermal energy is given by 

Eth = | / PdV (A12) 

= 6^a^poc^(^e')^^^_ (A13) 
= 296poc^a^. (A14) 
The magnetic field is expressed as 5 = 27raG^/^I] in our models, where the central surface density is given by 

/'''"max 

S = 2 / pdr (A15) 
Jo 

= 2poa / e-'^d^^ (A16) 



= 5.38poa. (A17) 

The magnetic energy is therefore given by 

^mag=/ —dV (A18) 
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= fpocyaYn... (/'""^^ e-dC j (A19) 
= 4061a^pocW- (A20) 



The kinetic energy is evaluated as 



Ekh 



pv^dV. 



(A21) 



The initial velocity field is generated using a random number. When ten velocity fields are generated by changing 
the the seed of the random number, we obtained £^kin in the range of (85.1 — 419) poc^a^ and with an average of 
220pocla^M^ . For the velocity field used in this paper, the kinetic energy is evaluated as ^kin = l^dpoc^a^M^ . 

Introducing the parameters of the density enhancement factor / (see equation [2]), we obtain scaling lows of 
Eth ^ ^kin (X f/^M^ ^grav ^ f^^ ^mag ^ f^^^^. which yield energy ratios of ^th/|^grav| = 0.836/-\ 

^kin/l^gravl = 0.394f-^M^ and ^mag/l^gravl = H.ba^. 

B. SURFACE-TO-VOLUME RATIO AND AXIS RATIOS 

In order to estimate the complexity of the density structure, the normalized surface-to-volume ratio is calculated. 
We calculate the surface and volume of a region ft{p) where the density is larger than a given threshold p and the 
point of maximum density is included. Based on Minkowski functional analysis, the volume V{p) and surface S{p) are 
calculated by 

V{p)=NoAV, 
S{p) = {6No-2Ni)AS, 

_ J 1 for a cell in ^{p) 
^ij,k ~ 1 otherwise ' 

where the subscripts i^j^k denote the cell numbers in the x, and z-directions, respectively. The symbols AV and 
AS* denote the volume and surface of a cell, respectively. 

The surface-to- volume ratio S{p)/V{p)'^^^ has a large value when the region ^{p) has a complex shape. The ratio 
S{p) /V{pY^^ has a minimum value when the region Vt{p) is spherical. We normalize the surface-to- volume ratio so 
that its value is unity for a sphere. Because the computational cell is cubic, the surface area of a sphere is given 
approximately by 67rr^ instead of 47rr^, where r denotes the radius of the sphere. The volume of the sphere is 
approximated by 47rr^/3. The normalized surface-to- volume ratio is therefore given by 



(Bl) 
(B2) 

(B3) 
(B4) 

(B5) 



R(P) = 32/36,1/3 v{p)V^ = O-l^^TW^- (^^^ 
We also measure the lengths of the principal axes of the region r^(p), calculating the moment of the coordinates 



where 7"i,j,/c = (^i,j,fe,i, n,j,/c,2, n,j,/c,3) = {xij^k^yi,j,k^ ^ij^k) represents the coordinates of the cell. The coordinates of 
the baricenter Tg are estimated by 

" Wp) ^ ^^d^kVij^kAV. (B8) 



The three principal axes are defined by the square roots of three eigenvalues of , arranged in ascending order 

(ai < a2 < as). The axis ratios are estimated as a2/ai and aa/ai. 

When the region Vt overlaps the boundary of the computational domain, care must be taken regarding the periodic 
boundary conditions. Moreover, the axis ratios are obtained only when the one-dimensional length of the region Vt is 
smaller than the width of the computational domain. 
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C. ANALYSIS OF DENSE REGION 

We calculate the mean values of the dense region of the cloud cores, p, J, f2, and the direction of the minor axis. 
The mean density of the dense region is defined by p = M/y, where 

M= f p{r)dV, (CI) 

Jp>0.1p 

max 

and 

V = [ dV. (C2) 

The integral /^^^ ^ dV denotes a volume integral over the region with p > O.lpmax- The mean magnetic field in 
the dense region is defined by 

B=^ [ B{r)dV. (C3) 

max 

The angular momentum in the dense region is defined by 

J 

/p>0.1p, 

where Vg and Vg denote the coordinates and velocity of the baricenter, estimated respectively as 



f (r - Vg) X [v{r) - Vg] p{r)dV, (C4) 

^P>0.1Pmax 

ty of the baricenter, estimated respectively as 
rg = ^ [ rp{r)dV, (C5) 

^P>0.1pmax 



and 

vg = ^ [ v{r)p{r)dV. (C6) 

The mean angular velocity in the dense region is defined by 

Ci = r^j, (C7) 

where / indicates the moment of inertia, of which the components are estimated as 

hj= [ [ir-rs)%-{ri-r,,i){rj-rgj)]p{r)dV, (C8) 

max 

and the subscripts z, j represent coordinate labels, i.e., x = ri, ^ = r2, z = rs, and Sij denotes the Kronecker delta. 
The axes of the dense region are derived from the moment of the coordinates, 

Kij = [ {Vi- rg,i){rj - rgj)p{r)dV. (C9) 

max 



The eigenvectors of Ki^j indicate the orientation of the principal axes. 
The velocity dispersion is given by 



(A^) = (i- / [v{r)-vgfp{r)dv\ 



1/2 



max 

The velocity dispersion for the radial velocity is also estimated as 

-,1/2 



1 

M 



Vr{rYp{r)dV 

p>0.1pmax 



(CIO) 



(Cll) 



where Vr{r) denotes the radial component of the relative velocity, v{r) — Vg, with respective to the position of the 



baricenter, Vg. 
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